import math

_longitude = 0
_latitude = 0
_utc_offset = 8
_day_of_year = 0
_delta = 0
_LSTM = 0
_TC = 0


def set_parameters(longitude=0, latitude=0, utc_offset=8, day_of_year=1):
    global _longitude, _latitude, _utc_offset, _day_of_year, _delta, _LSTM, _TC
    _longitude = longitude
    _latitude = latitude
    _utc_offset = utc_offset
    _day_of_year = day_of_year
    _delta = 23.45 * math.sin(math.radians(360 * (284 + _day_of_year) / 365))
    _LSTM = 15 * _utc_offset  # degrees
    B = math.radians(360 * (_day_of_year - 81) / 365)  # rad
    EOT = 9.87 * math.sin(2 * B) - 7.53 * math.cos(B) - 1.5 * math.sin(B)  # minutes
    _TC = 4 * (_longitude - _LSTM) + EOT  # minutes


def solar_hour_angle(local_time):
    LST = local_time + _TC / 60  # hours
    HRA = 15 * (LST - 12)  # degrees
    return HRA


def solar_elevation(local_time):
    HRA = solar_hour_angle(local_time)
    latitude_rad = math.radians(_latitude)
    delta_rad = math.radians(_delta)
    HRA_rad = math.radians(HRA)
    return math.degrees(
        math.asin(
            math.sin(latitude_rad) * math.sin(delta_rad)
            + math.cos(latitude_rad) * math.cos(delta_rad) * math.cos(HRA_rad)
        )
    )


def get_sunrise_sunset_time():
    latitude_rad = math.radians(_latitude)
    delta_rad = math.radians(_delta)
    cos_HRA = -math.tan(latitude_rad) * math.tan(delta_rad)
    HRA = math.acos(cos_HRA)
    HRA = math.degrees(HRA)
    sunrise = 12 - HRA / 15 - _TC / 60
    sunset = 12 + HRA / 15 - _TC / 60
    return sunrise, sunset


def get_noon_time():
    return 12 - _TC / 60


if __name__ == "__main__":
    import datetime as dt
    import numpy as np
    import matplotlib.pyplot as plt

    set_parameters(
        longitude=97.45979,
        latitude=28.65630,
        utc_offset=8,
        day_of_year=dt.datetime.now().timetuple().tm_yday,
    )
    sunrise, sunset = get_sunrise_sunset_time()
    sunrise_h = math.floor(sunrise)
    sunrise_m = (sunrise - sunrise_h) * 60
    sunset_h = math.floor(sunset)
    sunset_m = (sunset - sunset_h) * 60
    print(f"Sunrise: {sunrise_h}:{sunrise_m:.0f}, Sunset: {sunset_h}:{sunset_m:.0f}")
    noon = get_noon_time()
    noon_h = math.floor(noon)
    noon_m = (noon - noon_h) * 60
    print(f"Noon: {noon_h}:{noon_m:.0f}")
    time = np.linspace(0, 24, 1000)
    elevation = [solar_elevation(t) for t in time]
    plt.grid()
    plt.plot(time, elevation)
    plt.xlabel("Time (hour)")
    plt.ylabel("Solar Elevation (degree)")
    plt.title("Solar Elevation")
    plt.show()
